print("Don't run me")
# Import subset from NCCS BMF Data on nonprofits in Syr MSA
source_data("https://github.com/lecy/analyzing-nonprofit-service-areas/blob/master/DATA/syr_orgs.Rda?raw=true")
# Subset - Extract poverty orgs
pov_orgs <- syr_orgs[which ( syr_orgs$NTEECC == "B61" |
syr_orgs$NTEECC == "B82" |
syr_orgs$NTEECC == "B92" |
syr_orgs$nteeFinal1 == "E" |
syr_orgs$nteeFinal1 == "F" |
syr_orgs$nteeFinal1 == "I" |
syr_orgs$NTEECC == "K30" |
syr_orgs$NTEECC == "K31" |
syr_orgs$NTEECC == "K34" |
syr_orgs$NTEECC == "K35" |
syr_orgs$NTEECC == "K36" |
syr_orgs$nteeFinal1 == "L" |
syr_orgs$nteeFinal1 == "O" |
syr_orgs$nteeFinal1 == "P" |
syr_orgs$nteeFinal1 == "R" |
syr_orgs$nteeFinal1 == "S" |
syr_orgs$nteeFinal1 == "T"), ]print("Don't run me")
# download specific ggmap, register Google key
devtools::install_github("dkahle/ggmap")
register_google(key = 'AIzaSyAMrECdXqpOt443ms6nzh118uUsqC2lE6M',
account_type = "premium", day_limit = 15000)
# subset orgs for gecoding address information
pov_orgs_gps <- subset(pov_orgs, select=c(EIN:INCOME))
pov_orgs_gps$whole_address <- do.call(paste, c(pov_orgs_gps[
c("ADDRESS", "CITY","STATE",
"ZIP5")], sep = ", "))
pov_orgs_gps <- subset(pov_orgs_gps, select=c(whole_address, EIN:INCOME))
for (i in 1:nrow(pov_orgs_gps)) {
latlon = geocode(pov_orgs_gps[i,1])
pov_orgs_gps$lon[i] = as.numeric(latlon[1])
pov_orgs_gps$lat[i] = as.numeric(latlon[2])
}
# Make sure you find a way to then take any PO Boxes and place them in the
# center of centroids, and only the ones that are NAs.Download & clean Census ACS 2011-2015 5-yr data
print("Don't run me")
# Census API
# install.packages("devtools")
# devtools::install_github("hrecht/censusapi")
censuskey <- "f8064a17832b439e690be95ab0e0ef9a78617584"
var.list <- c("NAME",
"B25077_001E", "B25003_001E", "B25003_002E", "B25003_003E",
"B25002_001E", "B25002_002E", "B25002_003E", "B01003_001E",
"B19013_001E", "B17020_002E", "B12006_001E", "B23025_005E", "B15003_017E", "B15003_018E",
"B03001_001E", "B03002_003E", "B03002_004E", "B03002_005E",
"B03002_006E", "B03002_007E", "B03002_008E", "B03002_009E",
"B03001_002E", "B03001_003E",
"GEOID" )
acs_2015 <- getCensus( name="acs5",
vintage = 2015,
key = censuskey,
vars= var.list,
region ="tract:*",
regionin ="state:36"
)
# subset to Syracuse
# sum(acs_2015$county == "067")
# sum(acs_2015$county == "075")
# sum(acs_2015$county == "053")
acs_2015_syr <- acs_2015[which (acs_2015$county == "067" |
acs_2015$county == "075" |
acs_2015$county == "053"), ]print("Don't run me")
# Race variables are single race, non-hispanic. For example, "white" is really "white non-hispanic." This is true for every race category listed except, obviously, hispanic.
rename( acs_2015_syr,
c("B25077_001E" = "mdn_hous_val",
"B25003_001E" = "tenure_tot",
"B25003_002E" = "tenure_own",
"B25003_003E" = "tenure_rent",
"B25002_001E" = "tot_occup",
"B25002_002E" = "occupied",
"B25002_003E" = "vacant",
"B01003_001E" = "tot_pop",
"B19013_001E" = "mhh_income",
"B17020_002E" = "poverty",
"B12006_001E" = "labor_part",
"B23025_005E" = "unemploy",
"B15003_017E" = "high_sch",
"B15003_018E" = "ged",
"B03001_001E" = "tot_race",
"B03002_003E" = "white",
"B03002_004E" = "black",
"B03002_005E" = "am_ind",
"B03002_006E" = "asian",
"B03002_007E" = "islander",
"B03002_008E" = "other",
"B03002_009E" = "mixed",
"B03001_002E" = "not_hispanic",
"B03001_003E" = "hispanic"
))
acs_2015_syr <- rename.vars(acs_2015_syr,
c("B25077_001E", "B25003_001E", "B25003_002E", "B25003_003E",
"B25002_001E", "B25002_002E", "B25002_003E", "B01003_001E",
"B19013_001E", "B17020_002E", "B12006_001E", "B23025_005E", "B15003_017E", "B15003_018E",
"B03001_001E", "B03002_003E", "B03002_004E", "B03002_005E",
"B03002_006E", "B03002_007E", "B03002_008E", "B03002_009E",
"B03001_002E", "B03001_003E"),
c("mdn_hous_val", "tenure_tot", "tenure_own", "tenure_rent",
"tot_occup", "occupied", "vacant", "tot_pop",
"mhh_income", "poverty", "labor_part", "unemploy",
"high_sch", "ged",
"tot_race", "white", "black", "am_ind",
"asian", "islander", "other", "mixed",
"not_hispanic", "hispanic"))
# Create census race proportion variable
acs_2015_syr$porp_white <- (acs_2015_syr$white) / (acs_2015_syr$tot_race)
acs_2015_syr$porp_m <- abs((acs_2015_syr$porp_white) - 1)
names( acs_2015_syr )file.url <- "https://raw.githubusercontent.com/lecy/analyzing-nonprofit-service-areas/master/SHAPEFILES/syr_merged.geojson"
syr <- geojson_read( file.url, method="local", what="sp")
spTransform(syr, CRS("+proj=longlat +datum=WGS84"))## class : SpatialPolygonsDataFrame
## features : 186
## extent : -76.64471, -75.23954, 42.72384, 43.70702 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## variables : 68
## names : STATEFP.1, COUNTYFP.1, TRACTCE.1, GEOID.1, NAME.1, NAMELSAD.1, MTFCC.1, FUNCSTAT.1, ALAND.1, AWATER.1, INTPTLAT.1, INTPTLON.1, STATEFP.2, COUNTYFP.2, TRACTCE.2, ...
## min values : 36, 053, 030101, 36053030101, 301.01, Census Tract 301.01, G5020, S, 102677781, 0, +42.7960459, -075.3432945, 36, 067, 000100, ...
## max values : 36, 053, 031100, 36053031100, 311, Census Tract 311, G5020, S, 8864612, 93115, +43.1261230, -075.8928470, 36, 067, 940000, ...
file.url <- "https://raw.githubusercontent.com/lecy/analyzing-nonprofit-service-areas/master/SHAPEFILES/syr_merged_cen.geojson"
centroids <- geojson_read( file.url, method="local", what="sp")
spTransform(centroids, CRS("+proj=longlat +datum=WGS84"))## class : SpatialPointsDataFrame
## features : 186
## extent : -76.55702, -75.34358, 42.78044, 43.63666 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## variables : 1
## names : dummy
## min values : 0
## max values : 0
file.url <- "https://raw.githubusercontent.com/lecy/analyzing-nonprofit-service-areas/master/SHAPEFILES/pov_orgs.geojson"
pov_orgs <- geojson_read( file.url, method="local", what="sp")
spTransform(pov_orgs, CRS("+proj=longlat +datum=WGS84"))## class : SpatialPointsDataFrame
## features : 1095
## extent : -76.58932, -75.27822, 42.7584, 43.58445 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## variables : 21
## names : whole_address, EIN, FIPS, NTEECC, FILER, ZFILER, NAME, ADDRESS, CITY, STATE, ZIP5, GEN, SUBSECCD, RULEDATE, FNDNCD, ...
## min values : 1 ARKIE ALBANESE AVE, MANLIUS, NY, 13104, 010587946, 36053, B61, N, N, 0SWEGO COUNTY DEPUTIES ASSOCIATION, 1 ARKIE ALBANESE AVE, BALDWINSVILLE, NY, 13027, 0000, 02, 000000, 00, ...
## max values : SUNY COLLEGE OF FORESTRY, SYRACUSE, NY, 13210, 952505709, 36075, T99, Y, Y, ZONTA INTERNATIONAL SYRACUSE, SUNY COLLEGE OF FORESTRY, WEST EDMESTON, NY, 13485, 9434, 92, 201607, 21, ...
par( mar=c(0,0,0,0) ) # drop plot marginsplot( syr )
plot( centroids, col="black", pch = 20, add = T)plot( syr )
plot( centroids, col="black", pch = 20, add = T )
plot( pov_orgs, col="red", lwd = 2, add = T )names(syr)## [1] "STATEFP.1" "COUNTYFP.1" "TRACTCE.1" "GEOID.1"
## [5] "NAME.1" "NAMELSAD.1" "MTFCC.1" "FUNCSTAT.1"
## [9] "ALAND.1" "AWATER.1" "INTPTLAT.1" "INTPTLON.1"
## [13] "STATEFP.2" "COUNTYFP.2" "TRACTCE.2" "GEOID.2"
## [17] "NAME.2" "NAMELSAD.2" "MTFCC.2" "FUNCSTAT.2"
## [21] "ALAND.2" "AWATER.2" "INTPTLAT.2" "INTPTLON.2"
## [25] "STATEFP" "COUNTYFP" "TRACTCE" "GEOID"
## [29] "NAME" "NAMELSAD" "MTFCC" "FUNCSTAT"
## [33] "ALAND" "AWATER" "INTPTLAT" "INTPTLON"
## [37] "GEOID_full" "NAME.3" "GEOID.3" "state"
## [41] "county" "tract" "mdn_hous_val" "tenure_tot"
## [45] "tenure_own" "tenure_rent" "tot_occup" "occupied"
## [49] "vacant" "tot_pop" "mhh_income" "poverty"
## [53] "labor_part" "unemploy" "high_sch" "ged"
## [57] "tot_race" "white" "black" "am_ind"
## [61] "asian" "islander" "other" "mixed"
## [65] "not_hispanic" "hispanic" "porp_white" "porp_m"
map <- leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY")
map %>% addProviderTiles(providers$CartoDB)syr_map <- leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
addPolygons(data = syr, fillOpacity = 0.3) %>%
addCircles(data = pov_orgs) %>%
addCircles(data = centroids)
syr_map %>% addProviderTiles(providers$CartoDB) syr_map <- leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
addPolygons(data = syr, fillOpacity = 0.3, weight = 2) %>%
addCircles(data = pov_orgs, color= "red", opacity = 1.0) %>%
addCircles(data = centroids, color= "black", opacity = 1.0)
syr_map %>% addProviderTiles(providers$CartoDB) # change variables to be numeric
syr$porp_white <- as.numeric(as.character(syr$porp_white))
syr$porp_m <- as.numeric(as.character(syr$porp_m))
syr$mhh_income <- as.numeric(as.character(syr$mhh_income))
class(syr$porp_m)## [1] "numeric"
class(syr$porp_white)## [1] "numeric"
class(syr$mhh_income)## [1] "numeric"
# By race
pal <- colorNumeric(
palette = "Blues",
domain = syr$porp_m )
syr_map <- leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
addPolygons(data = syr, fillOpacity = 0.3, weight = 2,
color = ~pal(syr$porp_m)) %>%
addCircles(data = pov_orgs, color= "red", opacity = 1.0) %>%
addCircles(data = centroids, color= "black", opacity = 1.0)
syr_map %>% addProviderTiles(providers$CartoDB)# By Income level
pal2 <- colorNumeric(
palette = "Blues",
domain = syr$mhh_income )
syr_map2 <- leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
addPolygons(data = syr, fillOpacity = 0.3, weight = 2,
color = ~pal2(syr$mhh_income)) %>%
addCircles(data = pov_orgs, color= "red", opacity = 1.0) %>%
addCircles(data = centroids, color= "black", opacity = 1.0)
syr_map2 %>% addProviderTiles(providers$CartoDB)print("Don't run me")
# Create buffer, calculate number of nonprofit points around centroids
centroids_buff <- gBuffer(spgeom = centroids, width = .024)
spTransform(centroids_buff, CRS("+proj=longlat +datum=WGS84"))
CRS.new<-CRS("+proj=longlat +datum=WGS84")
proj4string(centroids_buff) <- CRS.new
proj4string(pov_orgs) <- CRS.new
over <- over( pov_orgs, centroids_buff)
# syr$locate <- over$GEOID_full
# syr$count <- rep.int(1, 186)
# other option?
# syr$test <- 1:nrow(over_clean)
pal2 <- colorNumeric(
palette = "Blues",
domain = syr$mhh_income )
syr_map2 <- leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
addPolygons(data = syr, fillOpacity = 0.3, weight = 2,
color = ~pal2(syr$mhh_income)) %>%
addCircles(data = pov_orgs, color= "red", opacity = 1.0) %>%
addCircles(data = centroids, color= "black", opacity = 1.0) %>%
addPolygons(data = centroids_buff, fillOpacity = 0.3, weight = 2,)
syr_map2 %>% addProviderTiles(providers$CartoDB)
# Not quite sure what is needed. How do I add up values for every centroid from the overlay?
# Similar to above, but combines over and removing NAs into one step
# pov_orgs_t <- pov_orgs
# pov_orgs_t$test <- pov_orgs[centroids_buff,]# Poportion of census tract that are high/low diversity
# number of nonprofits by type in high/low racially diverse areas
# number of nonprofits by type in low/high income census tracts
# Corr analysis between nonprofit access and census tract attributes
# Number of nonprofits location by census tract income or race
# Race
# Income